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Abstract. We present a novel computational framework for diffusive-reactive systems that sat- 
isfies the non-negative constraint and maximum principles on general computational grids. The 
governing equations for the concentration of reactants and product are written in terms of tensorial 
diffusion-reaction equations. We restrict our studies to fast irreversible bimolecular reactions. If 
one assumes that the reaction is diffusion-limited and all chemical species have the same diffusion 
coefficient, one can employ a linear transformation to rewrite the governing equations in terms of 
invariants, which are unaffected by the reaction. This results in two uncoupled tensorial diffu- 
sion equations in terms of these invariants, which are solved using a novel non-negative solver for 
tensorial diffusion-type equations. The concentrations of the reactants and the product are then 
calculated from invariants using algebraic manipulations. The novel aspect of the proposed com- 
putational framework is that it will always produce physically meaningful non-negative values for 
the concentrations of all chemical species. Several representative numerical examples are presented 
to illustrate the robustness, convergence, and the numerical performance of the proposed computa- 
tional framework. We will also compare the proposed framework with other popular formulations. 
In particular, we will show that the Galerkin formulation (which is the standard single-field for- 
mulation) does not produce reliable solutions, and the reason can be attributed to the fact that 
the single-field formulation does not guarantee non-negative solutions. We will also show that the 
clipping procedure (which produces non-negative solutions but is considered as a variational crime) 
does not give accurate results when compared with the proposed computational framework. 



1. INTRODUCTION AND MOTIVATION 

Mixing of chemical species across plume boundaries has a major influence on the fate of reactive 
pollutants in subsurface flows. In many practical cases, the intrinsic rate of reaction is fast compared 
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with other relevant time scales and hence the reaction may be assumed instantaneous (e.g., see 
references [HE]). Mixing is commonly modeled as an anisotropic Fickian diffusion process, with the 
effective diffusion coefficient aligned with the flow velocity (termed the longitudinal hydrodynamic 
dispersion coefficient) being much larger than the transverse components. Moreover, the small- 
scale spatial variability of permeability in real aquifers leads to highly heterogeneous velocity fields, 
which means that the principal directions of the diffusion tensor will vary in space and will not 
be aligned with the numerical grid. It has also been demonstrated that accurate treatment of 
dispersive mixing is crucial for computing the large-time spatial extent of the contaminant plume, 
which is a very important measure of the contaminant risk and thus of great interest to government 
regulators and the general public (e.g., see Reference |3j). It is well-known that heterogeneity and 
anisotropy lead to irregular plume boundaries, which enhance mixing-controlled reactions through 
increasing the interfacial area of the plume. It is, therefore, crucial to capture heterogeneity and 
anisotropy in order to properly model reactive transport in hydrogeological systems. 

One can capture heterogeneity computationally through adequate mesh refinement or by em- 
ploying multiscale methods (for example, see references [4], \E\, [6J, \7\). Although simulating flow 
and transport in highly heterogeneous porous media is still an active area of research, we shall not 
address this aspect in this paper. We shall employ meshes that are fine enough to be able to resolve 
heterogeneity. This paper focuses on resolving anisotropy to be able to accurately predict the fate 
of chemical species. We shall model the spatial and temporal variation of chemical species through 
diffusion-reaction equations. 

Diffusive-reactive equations arise naturally in modeling biological [81 [9], chemical [10L lllj . 
and physical [12 systems. The areas of application of diffusion-reaction systems range from con- 
taminant transport |13j . semiconductors [14] , combustion theory |15j to biological population 
dynamics [16j . It is well-known that these type of equations can exhibit complex and well-ordered 
structures/patterns |17|, 112] . A lot of effort has also been put into the mathematical analysis 
of diffusive-reactive systems. Many qualitative mathematical properties (e.g., comparison prin- 
ciples, maximum principles, the non-negative constraint) have been obtained for these systems 
|18|, I19|, [20J . A detailed discussion on mathematical aspects of diffusive-reactive systems is beyond 
the scope of this paper and is not central to the present study. 

1.1. Fast bimolecular reactions. The main aim of this paper is to present a robust compu- 
tation framework to obtain physically meaningful numerical solutions for diffusive-reactive systems. 
We shall restrict to bimolecular fast reactions for which the rate of reaction is controlled by diffusion 
of the two chemical species. Fast irreversible bimolecular diffusion-reaction equations come under 
the class of diffusion-controlled reactions |21l |22j . Examples of such systems include ionic reactions 
occurring in aqueous solutions such as acid-base reactions |23L [24 j , polymer chain growth kinetics 
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catalytic reactions [26J and enzymatic reactions |27J. In such systems the formation of the 
product is much faster than the diffusion of the reactants. In this paper we are only concerned about 
homogeneous reactions (i.e., the reactants are all in the same phase). Modeling diffusive-reactive 
systems in which reactions are heterogeneous (i.e., the reactants are in different phases) by incorpo- 
rating the surface effects of the reactants (for example, see references j28L 129), 130], 131) 132) 133), 134] ) 
is beyond the scope of the present paper but will be considered in our future works. 

Some main challenges in solving a system of diffusion-reaction equations are as follows: 

(a) Anisotropy: Developing robust computational frameworks for highly anisotropic diffusive- 
reactive systems is certainly gaining prominence. However, caution needs to be exercised in 
selecting numerical formulations to avoid negative values for the concentration of the chemical 
species. Many popular numerical schemes such as the standard single-field formulation [35J, 
the lowest order Raviart-Thomas formulation [36], and the variational multiscale mixed for- 
mulation |37) 138] give unphysical negative values for the concentration even for pure diffusion 
equations |39] , Furthermore, mesh refinement (either /i-refinement [40j or p-refinement |41j) 
will not overcome the problem of negative values for the concentration. Keeping in mind about 
these concerns and developing a robust framework for diffusive-reactive systems is a challenging 
task. 

(b) Nonlinear ity: The equations governing for these systems are coupled and nonlinear (see equa- 



tions (2.2a)-(2.2d), which are presented in a subsequent section). The solutions to these 
diffusive-reactive systems can exhibit steep gradients [42] , and a robust numerical solution 
procedure should be able to handle such features. 

(c) Scale dependence: These type of systems can exhibit multiple spatial and temporal scales. For 
example, the diffusion process can be much slower than the chemical reactions. Therefore, the 
numerical techniques should be able to resolve these multiple scales. Compared to fast reactions, 
different solution strategies are required for moderate and slow reactions, which are typically 
easier to solve due to smaller gradients. Construction of adaptive numerical techniques that 
take the advantage of specific reaction kinetics and simultaneously satisfying the underlying 
mathematical properties is still in its infancy. 

(d) Bifurcations, physical instabilities and pattern formations: These systems are capable of ex- 
hibiting physical instabilities, and even chaos [18J. 

We shall overcome the first challenge by employing a novel non-negative solver that ensures 
physically meaningful non-negative values for the concentration of chemical species. We employ 
a transformation of variables, which will overcome the second and third challenges. We solve 
problems that do not exhibit physical instabilities and chaos. Although considerable progress 
has been madw in developing numerical solutions of diffusive-reactive systems |43L 44J, none of 



these studies addressed the difficulties in obtaining non-negative solutions especially under strong 
anisotropy, which is the main focus of this paper. The following systematic approach has been 
employed to achieve the desired goal. The reactions are assumed to be fast and bimolecular. 
The concentration of reactants and product are governed by tensorial diffusion-reaction equations, 
and by an appropriate stoichiometric relationship. The three coupled tensorial diffusion-reaction 
equations are rewritten in terms of invariants, which are unaffected by the reaction. This procedure 
results in two uncoupled tensorial diffusion equations in terms of the invariants. A robust non- 
negative solver is employed to solve these two uncoupled tensorial diffusion equations following the 
ideas outlined in references [39, 40J. The concentrations of the reactants and product are then 
calculated using simple algebraic relations given by the transformation that is employed. 

1.2. An outline of the paper. The remainder of the paper is organized as follows. In Section 
[2j we first present the governing equations for a diffusive-reactive system in a general setting. We 
later restrict the study to fast bimolecular reactions, and deduce the corresponding governing 
equations. In Section [3j we present a novel non-negative computational framework for solving 
diffusive-reactive systems. Section [4] provides a theoretical discussion on the optimization-based 
solver for enforcing maximum principles. Representative numerical results are presented in Section 
[5j and conclusions are drawn in Section [6} 

The standard symbolic notation is adopted in this paper. We shall denote scalars by lowercase 
English alphabet or lowercase Greek alphabet (e.g., concentration c and density p). We shall make 
a distinction between vectors in the continuum and finite element settings. Similarly, a distinction 
is made between second-order tensors in the continuum setting versus matrices in the context of the 
finite element method. The continuum vectors are denoted by lower case boldface normal letters, 
and the second-order tensors will be denoted using upper case boldface normal letters (e.g., vector 
x and second-order tensor D). In the finite element context, we shall denote the vectors using lower 
case boldface italic letters, and the matrices are denoted using upper case boldface italic letters (e.g., 
vector v and matrix K). It should be noted that repeated indices do not imply summation. (That 
is, Einstein's summation convention is not employed in this paper.) Other notational conventions 
adopted in this paper are introduced as needed. 

2. GOVERNING EQUATIONS: DIFFUSIVE-REACTIVE SYSTEMS 

We now present governing equations of a simple diffusive-reactive system in a rigid porous 
medium. Let Q C M. nd be an open bounded domain, where "nd" denotes the number of spatial 
dimensions. The boundary is denoted by dQ, which is assumed to be piecewise smooth. Mathe- 
matically, d£l := & — $7, where Cl is the set closure of fL A spatial point in Cl is denoted by x. 
The divergence and gradient operators with respect to x are, respectively, denoted by div[-] and 
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grad[-]. The unit outward normal to the boundary is denoted by n(x). Let t e]0,Z[ denote the 
time, where I is the length of time of interest. We shall write the governing equations for the 
individual chemical species and for the mixture on the whole using the mathematical framework 
provided by the theory of interacting continua [45J . 

We shall assume that there are n chemical species, which include both reactants and product(s). 
We shall denote the molar mass of the i-th species by Q . The mass concentration of the i-th species 
at the spatial position x and time t is denoted by (fi(x,t). The molar concentration of the i-th 
species is denoted by Cj(x, t) and is defined as follows: 

ipijx, t) 

Ci(x,t) := — (2.1) 

For each chemical species (i = 1, • • • , n), the boundary is divided into two complementary parts: 
and . is that part of the boundary on which concentration for the i-th species is prescribed, 
and is the part of the boundary on which the flux for the i-th chemical species is prescribed. 
For well-posedness and uniqueness of the solution, we require Tf n Tf = 0, rp U Ff = <9ft, and 
meas(rP) > 0. The governing equations for the fate of the i-th chemical species can be written as 
follows: 

- div[D(x) grad[^j]] = mj(x,i) + gi(x, ipi, •••,(/?„, t) inftx]0,Z[ (2.2a) 
p f (x,i) = ^ p (x,t) onrfx]0,X[ (2.2b) 
n(x).D(x)grad[y>i]=0?(x,t) onrfx]0,X[ (2.2c) 
yj i (x,t = 0) = v?(x) in ft (2. 2d) 

where D(x) is the diffusivity tensor, (p? (x, t) is the prescribed mass concentration on the boundary, 
i?f(x, t) is the prescribed mass concentration flux on the boundary, ^(x) is the (prescribed) initial 
condition, <7i(x, ipi, • • • , cp n , t) is the rate of production/depletion of the i-th species due to chem- 
ical reactions, and TOj(x,i) is the prescribed external volumetric mass supply for the i-th species. 



Equations (2.2a)-(2.2d) are, respectively, the balance of mass, the Dirichlet boundary condition, 
the Neumann boundary condition, and the initial condition for the i-th chemical species. The 
governing equation for the balance of mass for the mixture on the whole takes the following form: 

n 

5^ 5i (x, pi,--- ,<p n ,t) = VxGft, Vt G]0,J[ (2.3) 

i=l 



In Engineering and Applied Sciences, the above system of equations (2.2a)— (2. 2d) is commonly 
referred to as time-dependent or non-stationary diffusion-reaction equations. In the theory of 
partial differential equations, the above system of equations is referred to as coupled semilinear 
parabolic partial differential equations |18j . These types of equations are capable of exhibiting 
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many interesting and complex features like instabilities, bifurcations, and even chaos (for example, 
see Pao [18J). However, in this paper we shall focus only on three interrelated mathematical 
properties that these type of equations satisfy: maximum principles, comparison principles, and 
the non- negative constraint. 

For completeness, we shall document the assumptions in obtaining the governing equations 



(2.2a)-(2.2d) under the mathematical framework offered by the theory of interacting continua. 
The main assumptions are as follows: 

(a) The porous medium in which the reaction takes place is assumed to be rigid. 

(b) The flow aspects in the porous media are neglected. Hence, the advection velocity of each 
chemical species in the mixture is zero. 

(c) The chemical species do not take partial stresses. 

(d) The diffusion process is based on the Fickian model. 

(e) The diffusivity tensor D(x) is the same for all the chemical species, which is the case with many 
diffusive-reactive systems. (However, it should be noted that in some chemical and biological 
diffusive-reactive systems the diffusivity tensor D(x) will be different for different species.) 

One can obtain a hierarchy of more complicated models by relaxing/generalizing one or more of the 
aforementioned assumptions, which will be considered in our future works. We shall now restrict 
our study to bimolecular reactions. 

2.1. Bimolecular reactions. Consider two chemical species A and B that react irreversibly 
to give the product C according to the following stoichiometric relationship: 

n A A + n B B -> n c C (2.4) 

where tia, ng and nc are (positive) stoichiometric coefficients. The rate of the chemical reaction is 
denoted by r, which (in general) can be a nonlinear function of the concentrations of the chemical 
species involved in the reaction. For a bimolecular reaction, one can write the rate of the reaction 
at the spatial point x and time t as follows: 

r(x, ca(x, t), cb(x, t), cc(x, t),t) = 9 -^— = ^— = H — ^— (2.5) 

kaQa n B Q B n c Cc 

Remark 2.1. One can handle stoichiometry using a stoichiometric matrix (for example, see 
references |45L I46j ). This approach will be convenient and appropriate for complicated chemical re- 
actions that involve multiple stages, multiple chemical species, and intermediate complexes. Herein, 
we shall not take such an approach, as the reaction is relatively simpler. 
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For bimolecular reactions, the governing equations (2.2a)-(2.2d) can be rewritten as follows: 

dc A 



dt 
dc B 

dt 
dc c 

dt 



div[D(x)grad[c A ]] = f A {x,t)-n A r mttx]0,l[ 
div[D(x)grad[c B ]] = /s(x,i) -n B r in fix]0,X[ 
div[D(x) grad[cc]] = fc( x , t) + nc r inf2x]0,Z[ 



Ci(x,t) = c?(x,t) onrfx]0,X[ (i = A, 5, C) 

n(x) • D(x) gradfei] = fcffoi) onrfx]0,X[ (i = A, B, C) 

Cj(x, t = 0) = c°(x) in 17 (t = A, 5, C) 



(2.6a) 

(2.6b) 

(2.6c) 
(2.6d) 
(2.6e) 
(2.6f) 



where the prescribed molar supply and the prescribed molar initial condition for the i-th species 
are, respectively, defined as follows: 

mi(x,t) 



/i(x,i) 
c°(x) 



Ci 



(2.7) 
(2.8) 



and the prescribed molar concentration and the prescribed molar concentration flux for the i-th 
species on the boundary are, respectively, defined as follows: 



/i?(x,t) 







(2.9) 
(2.10) 



3. PROPOSED COMPUTATIONAL FRAMEWORK 

We now present a robust numerical framework to solve bimolecular diffusive-reactive systems 



given by equations (2.6a)-(2.6f ) that will produce physically meaningful non-negative solutions for 



the concentration of chemical species on general computational grids. As mentioned earlier, we 
shall restrict our studies to fast reactions. By fast reactions we mean that the rate of the chemical 
reaction is much faster than the rate of the diffusion process. 

3.1. A non-negative invariant set, and uncoupled equations. Using an appropriate 
linear transformation, one can obtain two independent invariants, which are unaffected by the 
chemical reaction. As we shall see later that these invariants will be governed by uncoupled linear 



diffusion equations. For the bimolecular diffusive-reactive system (2.6), one can construct three 



different sets of independent invariants. Herein, we shall employ the set of independent invariants 

that inherits the non-negative property, and the corresponding linear transformation can be written 
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as follows: 



c F := c A + ( — ) c c (3.1a) 
\n c j 

cg ■= Cb + ( — J c c (3.1b) 
\n C J 

Remark 3.1. Let us denote the other two other sets of independent invariants by (U,V) and 
(M,N), which can be devised as follows: 



cu ■= ca - — cb, cy •■= ca+ — Cc 
n B J \n c ) 

'n B \ ( n B \ 

CM ■= cb - \ — ca, c N := c B + — c c 
™Aj \ncj 



(3.2a) 
(3.2b) 



As one can see from the above expressions, cjj and cm need not be always non-negative, which is not 
unphysical as U and M do not refer to any chemical species. They are just convenient mathematical 
quantities, and the possibility of negative values for these quantities is just a consequence of their 
mathematical definitions. We shall not use the above invariants as they do not possess the non- 
negative property, and hence are not appropriate to be used in a non-negative solver for diffusion- 
type equations. 



In the remainder of this paper, we shall assume that 



1 B 



Tq, and shall denote them 



by r D for notational convenience. Similarly, we shall assume that 



rN 
1 B 



Tq, and shall 



denote them by T N . These assumptions will facilitate application of the linear transformation 
to the prescribed boundary conditions c?(x, t) and tij-(x, i) (i = A, B, C). Using straightforward 



manipulations on equations (2.6a)-(2.6f ) based on the non- negative invariant set (3.1 ) and appealing 



to the above assumptions, one can show that the governing equations for the invariant F can be 
written as follows: 



dcF 
~dt~ 



div[D(x)grad[c F ]] = /j?(x,i) in nx]0, J[ 



c F (x,t) = c p F {x,t) := c P A {x,t) + I c c (x,t) onr D x]0,I[ 

n(x) • D(x) grad[c F ] = h p F {x,t) := h p A (x,t) + h p (x,t) onr N x]0,I[ 

CF (x, t = 0) = 4(x) := c°t(x) + f^j c c (x) in £1 



(3.3a) 
(3.3b) 

(3.3c) 

(3.3d) 



Similarly, the governing equations for the invariant G take the following form: 

dc G 



dt 



div[D(x) grad[c G ]] = / G (x,t) infix]0,X[ 



CG(x,t) = (x,t) := 4 (x,t) + eg, (x,t) on T D x]0,X[ 

n(x) • D(x) grad[c G ] = (x, t) := h P B (x, i) + 
cg(M = 0) = c° G (x) := c° B (x) + ( M") c c(x ) in n 



— W(x,i) onr N x]0,X[ 
n c J 



(3.4a) 
(3.4b) 

(3.4c) 

(3.4d) 



In the above equations, the prescribed molar supplies for the invariants F and G in terms of the 
chemical species A, B, and C take the following form: 



f F (x,t) = f A ( x ,t) + 
/G(x ) t) = /fl(x,t) + 



n c 
n G 



/c(x,t) 
/c(x,t) 



(3.5a) 
(3.5b) 



3.2. Specializing on fast bimolecular reactions. Recall that by a fast reaction we mean 
that the time scale of the reaction is much faster than the time scales associated with the diffusion 
process of the chemical species. Hence, for bimolecular reactions that are fast, it is reasonable 
to assume that the species A and B cannot co-exist at any given instance of time at the same 



location x. Using equations (3.1a)-(3.1b), one can rewrite the whole problem for fast irreversible 
bimolecular fast reactions in terms of the conserved quantities cf(x, t) and cq(x, t) as follows: 



ca(x, t) = max 
cs(x,i) = 
cc(x,t) = 



c F (pc,t) 



c G (x,t), 

nBj 



tla) 



max 



cp(x,t) + 
) (c F (x,t) - c A (*,t)) 



— J co(x,t),0 

n B J 



(3.6a) 
(3.6b) 
(3.6c) 



It should be noted that the original coupled diffusive- reactive system (2.6) is nonlinear and the 
solution procedure to obtain the concentrations of A, B, and C is still nonlinear. This is because of 
the fact that max[-, •] is a nonlinear operator. However, the obvious advantage is that the resulting 
equations are much simpler to solve. 

Before we discuss about the non-negative constraint, we shall introduce relevant notation. Let 
Dt '■= O,x]0,I], and St '■= d£lx]0,I]. Let C m (K) denote the set of all m-times continuously 
differentiable functions on an open set K. The set of all continuous functions on the (set) closure 
of K is denoted by C{K). The set of all functions that are m-times continuously differentiable in 
xel! and Z-times continuously differentiable in t G]0,X] is denoted by C m,l (DT)- From the theory 



of partial differential equations we have the following non-negative lemma for (pure) transient 
diffusion- type equations. 

Lemma 3.2 (The non- negative constraint for transient diffusion equation). Let c(x, t) G C 2,1 (Dt)(~) 
C(Dt) such that 

dc 



dt 



div[D(x)grad[c]] > in Dt (3.7a) 
a n(x) • D(x)grad[c] + /3 c(x, t) > on S T (3.7b) 
c(x,t = 0)>0 in SI (3.7c) 

where «o > 0, (3q > ; and uq + (3q > on St- Then c(x, t) > in Dt unless it is identically zero. 
PROOF. A proof can be found in McOwen [47], Evans [19] or Pao [18] . □ 

Using the above lemma, we now prove that all the chemical species in a bimolecular fast 
reaction also satisfy the non- negative constraint, which will be the main result of the paper on the 
mathematical front. 

Theorem 3.3 (The non-negative constraint for transient bimolecular fast reactions) . Letci(x, t) G 
C 2 ^(D T ) n C(D T ) for i = A, B, C. If /i(x,i) > 0, c?(x,t) > 0, c°(x) > and hf(x,t) > 
{i = A, B, C) then for bimolecular fast reactions we have q(x, t) > Vx G & and t e]0,X]. 



Proof. Since Ci(x,t) G C 2,1 (Z>r) n C(D T ) (i = A, B, C), using equations (|3.1a|)~(|3Tlb|), it is 



evident that c^(x, t), cg(x, t) G C ' (Dt) D C(Dt)- Based on the hypothesis given in the theorem 



and using equations (3.3b)-(3.3d), (3.4b)-(3.4d) and (3.5a)-(3.5b), one can conclude that 



/ F (x,t) > 0, 4(x,t) > 0, h p F (x,t) > 0, c° F (x,t) > 
/ G (x,t) > 0, c P G (x,t) > 0, h p G (x,t) > 0, c G (x,t) > 



Using Lemma 3.2, one can conclude that 

cHx,i) > andc G (x,t) > Vx G Q, t e]0,Z] 



From the relations given in equations (3.6a) and (3.6b), it is evident that 



CA(x,i) > andc B (x,t) > VxeO, t G]0,X] 



From equations (3.6a) and (3.6c), one can conclude that 

cc(x,i) = 



Ill 



c G (x, i) 



(3.8a) 
(3.8b) 

(3.9) 

(3.10) 

(3.11) 
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In either case, we have 



c c (x,t) > Vx G O, t G]0,X] (3.12) 
This concludes the proof. □ 



Remark 3.4. One can show that even the steady-state solution of bimolecular fast reactions 
satisfy the non-negative constraint. A non-negative result similar to Lemma \3.S\ exists for steady- 
state scalar diffusion equation (for example, see Gilbarg and Trudinger [20\). Using such a result 
and following the arguments presented in Theorem \3.3\ one can show that the concentrations of the 
chemical species in a steady-state bimolecular fast reaction also satisfy the non-negative constraint. 



3.3. Numerical solution strategy. The proposed numerical solution strategy for simulat- 
ing fast irreversible bimolecular reactions boils down to solving two uncoupled tensorial diffusion 
equations. The procedure can be outlined as follows: 



Solve the uncoupled tensorial diffusion equations given by ( 3.3a )-( 3.3d) and (3.4a)-(3.4d) 
to obtain cf and cq- 



• Using the relations given in equations ( 3.6a )— ( 3.6c), compute the concentrations of reac- 
tants A and B, and product C. 

In the next subsection we shall outline a non-negative numerical solver for tensorial diffusion equa- 
tion, which will be used to obtain the invariants F and G. It should be noted that if one obtains 
non-negative values for the invariants, based on the second step in the aforementioned procedure, 
it is obvious that the concentrations for A, B and C will be non-negative even in the numerical 
setting. 



3.4. A non-negative solver. In references |48L I39L I40L I49j numerical methodologies have 
been proposed for diffusion-type equations based on optimization techniques. In this paper, we shall 
extend these methodologies to both steady-state and transient solutions for diffusion-controlled 
bimolecular- reactive systems. 

3.4.1. Steady-state analysis: The weak statements based on the standard Galerkin formulation 
to obtain the steady-state solutions for the invariants F and G can be written as follows: Find 
cf(x) G Vf and cg(x) G Vg such that we have 



B{q;c F ) = L F {q) Vg(x) G Q (3.13a) 

B{q;c G ) = L G {q) Vg(x) G Q (3.13b) 
ii 



where the bilinear form and the linear functionals are defined as follows: 

B(q; c) := / grad[g(x)] • D(x) grad[c(x)] dft (3.14a) 
Jn 

L F {q) ■= [ ?(x) / F (x) dQ + I q(x) (x) dT (3.14b) 
Jn Jr" 

L G {q) := f ? (x) / G (x) dft + / g(x) (x) dr (3.14c) 

In the above discussion, the following function spaces have been used: 

V F := (c(x) G H 1 ^) | c(x) = 4(x) on T D } (3.15a) 
V G := (c(x) G | c(x) = C p (x) on T D } (3.15b) 

Q := (g(x) G fT 1 ^) | g(x) = on T D } (3.15c) 

where H 1 ^) is a standard Sobolev space |50j . However, it should be noted that the Galerkin 
formulation does not satisfy the non-negative constraint and maximum principles for diffusion-type 
equations (see references [39, 40, 49J). These violations will be exacerbated in the case of diffusive- 
reactive systems (see Section[5]in this paper). We now modify the Galerkin formulation to meet the 
non- negative constraint and maximum principles for fast bimolecular diffusive-reactive systems. 



From Vainberg's theorem j5 1L [52], the weak form (3.13a) has an equivalent minimization 
problem, which can be written as follows: 

minimize -B(cf;cf) — Lf(cf) (3.16a) 
c F (x)eP F 2 

subject to cf n < cjr(x) < cf ax (3.16b) 

where c™ m and c™ ax are, respectively, the lower and upper bounds for the invariant F. These 
bounds arise due to the non-negative constraint and the maximum principle. For example, to 
enforce the non- negative constraint one can take c™ m = 0. 

After performing finite element discretization using low-order finite elements on the above 
minimization problem, a numerical methodology that satisfies the non-negative constraint and the 
maximum principle for the invariant F can be written as follows: 

minimize -(cf;Kcf) — (cp; fp) (3.17a) 

subject to cf in l <c F < cf ax l (3.17b) 

where (•; •) represents the standard inner-product on Euclidean spaces, 1 denotes a vector of ones 
of size ndofs x 1, and the symbol ^ represents the component-wise inequality for vectors. That is, 

a H b implies a» < 6j Vz (3.18) 
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The other notation employed in the above equations (3.17a)-(3.17b) is as follows: cp is the vector 
containing the nodal concentrations of the invariant F, f F is the nodal volumetric source vector, 
and K is the coefficient matrix, which will be symmetric and positive definite. The number of 
degrees-of- freedom in the nodal concentration vector is denoted by ndofs. Hence the size of each 
of the vectors cp and f F will be ndofs x 1, and the size of the matrix K will be ndofs x ndofs. 
A numerical methodology for the invariant G can be developed in a similar manner by replacing 
f F with the corresponding volumetric source vector f G , c™ m with c™ n and c™ ax with Cg ax . 

It should be noted that the coefficient matrix K (which is sometimes referred to "stiffness" 
matrix) for the invariant G is the same as that of the invariant F. This is due to the fact that 
the diffusivity tensor D(x) and the computational domain are the same for both the invariants. 



Since the coefficient matrix K is positive definite, the constrained optimization problem (3.17) 
belongs to convex quadratic programming and has a unique global minimizer [53] • The various 
steps in the proposed methodology to obtain steady-state solutions for bimolecular fast diffusive- 
reactive systems are summarized in Algorithm [TJ which will serve as a quick reference for computer 
implementation. 



Remark 3.5. It should be noted that the shape (or interpolation) functions for low- order finite 
are non-negative within the element. Therefore, enforcing non-negative constraints on the nodal 
concentrations will ensure non-negativeness within the element, and hence in the entire computa- 
tional domain. It should also be noted that the above methodology cannot be extended to high-order 
finite elements, as high-order interpolation functions can be negative within the finite element. For 
further details see references [4Ql 141] . 



Remark 3.6. A naive implementation of the numerical strategy proposed in subsection \3.S\ will 
lead to propagation of numerical noise from the invariant set to reactants and product. So care 
should be exercised while implementing the proposed numerical strategy to find the concentration 
of A, B, and C. As outlined in Algorithm^ first the numerical noise from invariants F and G 
should be removed so that product C does not get affected adversely. 

3.4.2. Transient analysis: We will first perform temporal discretization using the method of 
horizontal lines and then followed by spatial discretization. To this end, the time interval [0,1] is 
discretized into N non-overlapping sub-intervals: 

N 

[0,1] = lj[t n _i,t n ] (3-19) 

n=l 
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Algorithm 1 Numerical framework for steady-state fast bimolecular diffusive-reactive systems 



1: Input: Stoichiometric coefficients; volumetric sources; and boundary conditions for the chemical 
species A, B, and C 

2: Calculate corresponding volumetric sources and boundary conditions for the invariants 
3: Calculate c™ m , c™ ax , c™ m and cg ax based on maximum principles and the non-negative con- 
straint 

4: Call optimization-based steady-state diffusion solver to obtain c F : 

minimize -(c.f; Kc F ) - (c F ; f F ) 

c F eM. nd °f s ^ 

subject to cf n l ^c F < c£ ax l 



5: Call optimization-based steady-state diffusion solver to obtain c G : 

minimize -(c G ; Kc G ) - {c G ; f G ) 

c G £R nd °f s I 

subject to cg in l <c G < c£ ax l 



6: for i = 1, 2, ■ ■ ■ , Number of nodes in the mesh do 

7: if -e mach < c F (i) - c G (i) < +e mach then 

8: c A (i) = 0, c B (i) = 0, c c {i) = (^) c F (i) 

9: end if 

10: if c F (i) - c G (i) < -e mach then 

11: c A (i) = 0, c B (i) = - c F (i) + c G (i), c c {i) = c F (i) 

12: end if 

13: if c F (i) - c G {i) > +e mach then 

14: C A (i) = C F (i) - CG(i), C B (i) = 0, C C (i) = C G (i) 

15: end if 

16: end for 



where to = and tjv = I. For simplicity, we shall assume uniform time step, which will be denoted 
by At. That is, 

At := tn+i - t n Vn (3.20) 

However, it should be noted that a straightforward modification can handle non-uniform time 
steps. The following notation is employed to obtain time discretized version of quantities cf(x, t) 
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and 9cf q*'^ at t = t n : 



c^(x) := c F (x,t n ) 



dc F (yt,t n ) 



(3.21a) 
(3.21b) 



' dt 

Using the method of horizontal lines along with the backward Euler time stepping scheme, the 



governing equations ( 3.3a )-( 3.3d) at the time level i n +i can be rewritten as the following anisotropic 
diffusion with decay: 



1 \ (n+l)/ . r 
, c F '(x) - div 



Dfx) erad 



4 n+1) (x)H =/ F ( x ,t n+1 )+(^)4 n) (x; inO (3.22a) 



c^ +1) (x) = c^(x,t n+ i) := c^(x,i n+ i) + ( — ) c^(x,t n+ i) on T D 



nc J 



(3.22b) 



n(x) • D(x) grad 4" +1) (x)j = h* (x, t n+1 ) := (x, t n+1 ) + ( ~) h» (x, t n+1 ) on T N (3.22c) 



c F (x,t ) = 4(x) := 4(x) + ( ^ <&(x) in fi 



where f F (x,t n+ i) := f A {x,t n+1 ) + ( ^ ) /c(x, W:i 



(3.22d) 



Remark 3.7. From Section^ we have assumed that the volumetric sources, initial and bound- 
ary conditions for A, B, and C is non-negative. Correspondingly the the volumetric source /f(x, t n+ \) 
initial condition c F (x), the Dirichlet boundary condition c^(x, t n +\), and the Neumann boundary 
condition /i^(x, t n +\) are non-negative. As the decay coefficient in front of c" F +l \~x) is always 
greater than zero and the diffusivity tensor D(x) is assumed to be symmetric, bounded, and uni- 
formly elliptic, a maximum principle under these regularity assumptions for this type of anisotropic 
diffusion with decay exists |40j. 



The standard single-field formulation for equations (3.22a)-(3.22d) is given as follows: Find 

(3.23) 



c (n+i) ^ e gygh faofc we have 

B t (q ] c^ +l) )=L F {q) V 9 (x)GQ 
where the bilinear form and linear functional are, respectively, defined as 

B l (q; 4 n+1) ) := ^ J g(x) 4" +1) (x) dft + J grad[g(x)] • D(x) grad[4 n+1) (x)] dft (3.24a) 

L F (q) := jf g (x) (^/ F (x,t n+1 ) + ^4" } ( x )) ^ + ^ g(x) /£(x,t n+ i) dr (3.24b) 
The function space 7^ is defined as follows (for more rigorous definition see [19 s Section 7.1]): 

V F := {c F (;t) e H 1 ^) | c F (x,t) = c£(x,i) onT } (3.25) 
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but the function space Q denned by equation (3.15c) remains the same. In a similar manner, one 



can obtain a numerical methodology for transient solutions for the invariant G. 



An equivalent minimization formulation for the weak form (3.23) can be written as follows: 



minimize (c? +1) ; c£ +1) ) - L F (c { ; +1) ) (3.26a) 



> +1) (x)e^ 2 



subject to cf m < c^ +1) (x) < cf ax (3.26b) 

After finite element discretization by means of low-order finite elements, a non negative formulation 
based on the minimization problem can be devised as follows: 

minimize \ (c^-Kc^) - (c^hf^) - J- (c^;cW (3.27a) 
subject to cf n l < c ( p +1) < c£ ax l (3.27b) 

where c'p and f^ F +1 ^ are the nodal concentration and volumetric source vectors for invariant F 
at time level t n+ \. Correspondingly, the nodal concentration vector obtained at time level t n 
also acts as a volumetric source similar to f^p +1 \ As we have employed low-order finite elements 
for spatial discretization, we have c^) y and ^ 0. Hence at each time level, the net 

volumetric source vector + ^ 0. Note that the symmetric positive definite matrix K 

is different to that of K , which is the coefficient matrix in steady-state analysis. Since the coefficient 



matrix K is positive definite, the discrete minimization problem (3.27) belongs to convex quadratic 
programming, and from the optimization theory a unique global minimizer exists. Algorithm 
[2] outlines the steps involved in the implementation of the proposed methodology for transient 
diffusive-reactive systems. 

4. PROPERTIES OF THE PROPOSED COMPUTATIONAL FRAMEWORK: A 

THEORETICAL STUDY 

In a continuous setting, it is well-known from the theory of partial differential equations that the 
comparison principle, the maximum principle and the non-negative principle satisfy the following 
commutative diagram: 

uniqueness of the solution < comparison principle > maximum principle 



non-negative constraint 

For linear problems, one can also directly prove the uniqueness of the solution using the maximum 
principle (instead of the comparison principle). However, for semilinear and quasilinear elliptic 
and parabolic partial differential equations, the comparison principle is more convenient to prove 
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Algorithm 2 Numerical framework for transient fast bimolecular diffusive-reactive systems 



1: Input: Time step At; total time of interest I; stoichiometric coefficients; volumetric sources; 
initial and boundary conditions for the chemical species A, B, and C 

2: Calculate corresponding volumetric sources, initial and boundary conditions for the invariants, 
and form the vectors Cp and Cq 

3: Calculate c™ m , c™ ax , c™ n and cj^ ax based on maximum principles and the non-negative con- 
straint 

4: for n = 0,1,- ■ ■ ,N- 1 do 

5: Call optimization-based steady-state diffusion with decay solver to obtain c^" +1 ^: 
subject to cf n l < 4 n+1) < cf s 



Call optimization-based steady-state diffusion with decay solver to obtain c[? +1 ^: 

1 / ('n.-l-n ~ fn+1l\ / Cra-H"! .fn+1l\ 1 



c, 



minimize - K eg*' - cg'«>; - ±- <T "l c<?> 



7: 


for % = 1,2, • • ■ , 


8: 


11 ^mach C 


9: 


cj +1) (0 = 


10: 


end if 


11: 


if 4 n - 


12: 




13: 


end if 


14: 


if c? +1) (z) - 


15: 


cf +1) (0 = 


16: 


end if 


17: 


end for 



subject to cg in l ^ cg l+1) ^ cg ax ] 



«"(^j C G «<+ e mach then 

nc\ 



o, 4 n+1) « = o, 4 n+1) « = feU n+1) (*) 



=£) cg l+1) (i)<-e mach then 



cg +1 \i) > +e mach then 



18: end for 



directly the uniqueness of the solution. As mentioned earlier, the diffusive-reactive system given by 



equation (2.2) is a semilinear partial differential equation. This is the reason why we have indicated 
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in the above commutative diagram that the uniqueness of the solution is a direct consequence of 
the comparison principle. 

It is natural to ask which of these properties are inherited in discrete setting by the optimization- 
based solver for solving diffusion- type equations. Prior studies on the optimization-based method- 
ologies (e.g., references [481 I39L 1401 149] ) did not address whether such a commutative diagram 
is satisfied in the discrete setting. In particular, these studies did not discuss the discrete version 
of the comparison principle. The standard comparison principle takes the following form: Let 
C[u] := a(x)u(x) - div[D(x)grad[«]]. Let ui(x), u 2 (x) e C 2 (Q) n C°(Tl). If C[m] > C[u 2 ] in SI 
and ui(x) > u 2 (x) on d£l then iti(x) > it 2 (x) in fi. We shall say that a numerical formulation 
inherits the comparison principle (or possess a discrete comparison principle) if for any two finite 
dimensional vectors satisfying / 1 y / 2 implies the finite dimensional solutions under the numerical 
formulation satisfy c\ y c 2 . 

In discrete setting, the optimization-based methodology has the following properties: 

(i) The methodology is based on the finite element method. 

(ii) The methodology satisfies both the non-negative principle and the maximum principle on 
general computational grids and with no additional restrictions on the time step, 

(iii) The methodology does not satisfy the comparison principle. A counterexample is provided 
in Figure [T] However, if / 1 = rjf 2 with rj > 0, then c\ = nc 2 under the optimization-based 
solver. 

(iv) Although a discrete version of the comparison principle does not hold, a unique solution always 
exists under the optimization-based solver. The existence and the uniqueness of the solution 
in the discrete setting stems from the results in the convex programming with bounds on the 
variables. Of course, the proofs from the optimization theory do not require a discrete version 
of the comparison principle to hold but use alternate approaches from Matrix Algebra and 
Analysis |54j. 

(v) The solution procedure is nonlinear, as a quadratic programming optimization problem with 
inequalities is nonlinear. 

(vi) One can employ the consistent capacity matrix that stems from the variational structure of 
the underling weak formulation. If needed, one can also employ a lumped capacity matrix, 
which is considered to be a variational crime. 

In conclusion, the proposed optimization methodology does not possess a commutative diagram 
similar to the one given above. But it should be noted that the optimization-based method is the 
only known methodology that satisfies the non-negative principle and the maximum principle on 
general computational grids with no additional restrictions on the time step. Some other noteworthy 
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efforts on non-negative solvers for diffusion-type equations are by Le Potier [55J and Lipnikov et 
al. j56L 15 7j . which have the following attributes: 

(i) The methodology is based on the finite volume method. 

(ii) The methodology satisfies the non-negative principle on computational grids with triangles. 
The method does not satisfy the maximum principle. 

(iii) The methodology does not satisfy the comparison principle. 

(iv) The solution procedure is nonlinear, as the collocation points (i.e., the location of the unknown 
concentrations) in turn depend on the concentration. 

(v) One can employ only a lumped capacity matrix, which is consistent with the finite volume 
method. Also, one has to employ the backward Euler time stepping scheme, as any other time 
stepping scheme can result in an algorithmic failure or can produce negative concentrations. 



5. NUMERICAL PERFORMANCE OF THE PROPOSED COMPUTATIONAL 

FRAMEWORK 

In this section we shall illustrate the performance of the proposed computational framework 



for diffusive-reactive systems given by equations (2.6a)-(2.6f ). We shall also perform numerical 
/i-convergence using hierarchical meshes. We shall solve three representative problems: plume de- 
velopment from a boundary, plume formation due to continuous point source emitters, and decay 
of a chemically reacting slug. These types of problems have many practical applications. For ex- 
ample, these problems naturally arise from situations in which one may want to regulate/control 
an induced contaminant in an ecological system, which is often achieved by introducing a suitable 
chemical or biological species that reacts with the contaminant to form a harmless product. Some 
specific examples include chemical dispersion methods to control oil spills and airborne contam- 
inants, and usage of bioremediators for regulating industrial emissions into water bodies. In all 
the numerical simulations in this paper, the resulting convex quadratic programming problems are 
solved using the interior- point-convex quad prog algorithm available in MATLAB |58j. The toler- 
ance in the stopping criterion in solving convex quadratic programming is taken as 100e mac h, where 
e mach ~ 2.22 x 10 -16 is the machine precision for a 64-bit machine. 



5.1. Numerical /i-convergence study. We shall use the method of manufactured solutions 
to perform numerical /i-convergence. We shall restrict the current study to steady-state response. 
The computational domain is a rectangle of size L x x L y . The diffusivity tensor for this numerical 
study is taken as follows: 



D(x) = RD diagona iR 
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(5.1) 



where 



/ cos{6) -sin(0) \ , x 

R= W W 5.2 

\ sin(0) cos(0) / 

/ di \ 

Ddiagonal = (5.3) 

V d 2 / 

The assumed analytical solutions for cp(x) and cg(x) are taken as follows: 

Cf(x) = sin (S) sin (m9 (5 - 4a) 

cg(x) = cos (S) cos (i;) (5 - 4b) 

where < x < L x and < y < L y . Using the above expressions, the source terms, boundary 
conditions, and expressions for the concentrations of A, B, and C are in turn calculated. The 
corresponding volumetric sources for invariants F and G are given as follows: 

t , s vr 2 . / vrx \ . / vry \ / , /cos 2 ^) sin 2 (0)\ , /sin 2 (#) cos 2 (#) 

7r^ Z' vrx \ Z' Try \ 

- or r ( d i - d 2 ) sin((9) 003(61) cos — cos — (5.5a) 

t . , vr 2 /ttx\ / rry \ / /cos 2 (0) sin 2 (0)\ , /sin 2 ^) cos 2 ^) 

- 5T~T~ (* " 4)sin(»)c<M(»)sin (^) sin f^) (5.5b) 

A pictorial description of the boundary value problem is given in Figure [2] In this paper, we have 
taken the following values for the parameters to perform the steady-state numerical /i-convergence: 

L x = 2, L y = 1, 9 = 7t/3, di = 1000, d 2 = 1, n A = 2, n B = 3, n c = 1 (5.6) 

Figure [3] shows the typical computational meshes employed in this study. Figure [4] illustrates the 
convergence of the concentration of the invariants, and the convergence for the concentrations of the 
chemical species A, B, and C is shown in Figure [5] For invariants, we have shown the convergence 
in both L 2 norm and H 1 seminorm. However, for the chemical species A, B, and C, we have shown 
the convergence only in L 2 norm as max[-] operator is non-smooth. 

5.2. Plume development from boundary in a reaction tank. The test problem is pic- 
torially described in Figure |6j We restrict the present study to steady-state analysis. The stoichio- 
metric coefficients are taken as ua = 1, tib = 1 an d nc = 2. The diffusivity tensor is taken from 
the subsurface literature |13j . and can be written as follows: 

Dsubsurface(x) = a T ||v||I + ^jj^j-^V ® V (5.7) 
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where ® is the tensor product, v is velocity vector field, and ar and oll are, respectively, transverse 
and longitudinal diffusivity coefficients. It should be emphasized we have neglected the advection, 
and the velocity is used to define the diffusion tensor. The velocity field is defined through a 
multi-mode stream function, which takes the following form: 

3 



^(x,y) 



-y 



T, A 

k=l 



fcCOS 



PkTTX 



8,11 l^rJ 



(5.f 



The corresponding components of the velocity are given by 
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v x (x, y) 



(x,y) = + 



dip 
dy 

dip 
dx. 



1 + J> 



k=l 



k — COS 

Ly 



3 



A P k7r ■ 

/, A k-j— sin 
k=l 



L , 



( qkiry\ 



7T 



(5.9a) 



(5.9b) 



It is noteworthy that the above velocity vector field is solenoidal (i.e., div[v] = 0). Although 
realistic aquifers exhibit spatial variability on a hierarchy of scales, periodic or quasi-periodic models 



similar to the one outlined above (given by equations (5.7)-( |5.9| )) have often been used due to their 
simplicity [591160] . 

In this paper we shall take (L x ,L y ) = (2,1), (pi 5 p 2 ,P3) = (4,5,10), (91,92,93) = (1,5,10), 
(Ai, A2, A3) = (0.08,0.02,0.01), and (aL,ar) = (1,0.0001). The contours of the stream function, 
and the corresponding velocity vector field are shown in Figure [7| The computational domain is 
meshed using four-node quadrilateral elements. The numerical simulation is carried out on various 
structured meshes with XSeed and YSeed ranging from 97 to 259, where XSeed and YSeed are, 
respectively, the number of nodes along x-direction and y-direction. That is, a structured mesh 
with XSeed = YSeed = 259 has over 67,000 unknown nodal concentrations. 

The values for prescribed concentrations on the boundary are taken as <? A = 1 and c? B = 10. 
Figures [8 -12 compares the contours of the concentrations of the invariants, reactants and product 
under the Galerkin formulation, the clipping procedure and the proposed non-negative formulation. 



Figure 13 shows the regions in the domain that have violated the lower and upper bounds on the 



concentrations of the invariants and the product under the Galerkin formulation. To provide more 
quantitative inference on the violations, the variation of the concentration of the product along two 
sections given by y = 0.45 and x = 1.75 are plotted in Figure [14| 



Figures 15 and 16 show the variation of the total concentration, j c(x, y) dy, along x-direction for 
the invariants and the product. The total concentration for a given cross-section is a useful measure 
in the studies on reactive-transport. The results are shown for various mesh refinements, and the 
Galerkin formulation predicts unphysical negative values even for this important global measure. 

On the other hand, the proposed computational framework predicts physically meaningful for the 
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considered measure. Figures 17 and 18 show the contours of the Lagrange multipliers in solving 
quadratic programming problems to obtain the invariants F and G. The Lagrange multipliers arise 
due to the enforcement of the lower and upper bounds due to the non-negative constraint and the 
maximum principle in calculating the invariants. It should be noted that, based on the Karush- 
Kuhn- Tucker results in the theory of optimization, the Lagrange multipliers are non-negative 



5.3. Plume formation due to multiple stationary point sources. The problem is pic- 



torially described in Figure 19 The following rotated heterogeneous diffusivity tensor is employed: 

D(x) = RD (x)R T (5.10) 

where Dq(x) is 

D 0(X) = 9 9 ( 5 - 11 ) 

V-l-exy ey 2 +x 2 / 



and the rotation tensor is same as before (given by equation (5.2)) with angle 9 = ir/3. For this 



problem, the parameter e is taken as 0.001. The stoichiometry coefficients are taken to be ha = 1, 



riB = 1 and ric = 2. Figures 20-22 show the contours of the concentrations of the invariants 
and the product under the Galerkin formulation and the proposed computational framework. Fig- 
ure 



23 shows the variation of the integrated concentration of the product with respect to y (i.e., 
J c(x, y) dy) along x, and the Galerkin formulation dramatically predicts negative values even for 
this average. Table [T] shows that small violations in the non-negative constraint for invariants can 



lead to bigger violations in the non-negative constraint for the product C. Figure 24 compares the 
CPU time taken by the Galerkin formulation and the proposed computational framework. Even 
for a mesh with more than 40,000 nodes, the CPU time taken by the proposed framework is only 
1.07 times the CPU time taken by the Galerkin formulation. For this mesh, 50.51% of the nodes 
have violated the non-negative constraint for the product under the Galerkin formulation. On the 
other hand, the nodal concentrations of the product under the proposed computational framework 
are physically meaningful non-negative values. 

5.4. Diffusion and reaction of an initial slug. The test problem is pictorially described 
in Figure 25 The initial conditions within the reaction slug are assumed to be c A = 10, c B = 0, 



and c c = 0. In the other parts of the rectangular domain we have Ca = c b = c c = 0- ^he 
boundary conditions for the chemical species on the rectangular domain are assumed to be c A = 0, 
<? B = \ — e - *, and c c = 0. The problem is solved within the time interval t € [0, 1]. The diffusivity 
tensor within the domain is taken as D(x) = RD su b SU rfaceR T with 6 = ? . The stoichiometric 
coefficients are ua = 2, hb = 2, and ric = 1. Numerical simulations are performed using structured 
mesh with 101 nodes along each side of the domain using four-node quadrilateral elements. The 
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Table 1. Plume formulation due to multiple stationary point sources: This table shows the 
minimum concentration, maximum concentration, and % of nodes violated the non-negative 
constraint under the Galerkin formulation. The values clearly indicate that small violations 
in meeting the non-negative constraint for invariants can result in much bigger violation of 
the non-negative constraint for the product. 



Mesh 



Invariant F 



Min. cone. 
Max. 



cone. 



x 100 



nodes violated 



Invariant G 



Min. cone. 
Max. cone. 



x 100 



nodes violated 



21 x 21 
51 x 51 
101 x 101 
161 x 161 
201 x 201 



-1.50x10--' 

1.42x10° 
-1.38xlQ- 2 

2.37x10° 
-3.13xl0~ 2 

3.38x10° 
-4.45 xl0~ 2 

4.24x10° 
-4.98xlQ- 2 

4.71x10° 



x 100 
x 100 
x 100 
x 100 
x 100 



-0.11% 
-0.58% 
-0.93% 
-1.05% 
-1.06% 



8.62% 
24.14% 
25.30% 
25.77% 
27.29% 



-2.07X1Q- 3 



3.00x10" 


i 


-8.54x10 


-3 


4.99x10" 


1 


-1.31x10 


-2 


7.11x10" 


1 


-1.52x10 


-2 


8.96x10" 


1 


-1.58x10 


-2 



Mesh 



21 x 21 
51 x 51 
101 x 101 
161 x 161 
201 x 201 



Product C 



]f in - conc - x 100 nodes violated 

Max. cone 



-4.13X1Q- 3 



1.22x10- 


-l 


-2.75x10 


-3 


1.83x10- 


1 


-6.25x10 


-2 


2.26x10- 


1 


-8.90x10 


-2 


2.48x10" 


1 


-9.96x10 


-2 



x 100 = -3.38% 31.29% 

x 100 = -15.05% 53.94% 

x 100 = -27.70% 51.02% 

x 100 = -35.93% 54.94% 

x 100 = -38.90% 50.51% 



23.36% 
43.71% 
46.53% 
45.20% 
40. 



volumetric source of each species A, B, and C is equal to zero inside the rectangular domain. The 
time step is taken as At = 0.05 s. Figure [26| compares the concentrations of the product under the 



Galerkin formulation and the proposed computational framework. Figure 27 shows the variation 
of the concentration along x at y = 2.5 for t = 0.05 s, t = 0.1 s and t = 0.15 s. As evident from 
these figures, even for this transient problem, the proposed computational framework produced 
physically meaningful non-negative values for the concentration of the product at all time levels. 
Figure [28] shows the number of iterations taken by the quadratic programming solver at each time 
level to obtain concentrations of the invariants F and G. For this transient problem, the CPU times 
taken by the proposed computational framework and the Galerkin formulation are 1516.825 and 
1486.19. That is, the computational overhead of the proposed computational framework is 2% when 
compared with the Galerkin formulation. The elasped CPU time is measured using MATLAB's 
builtin tic-toe feature. 
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6. CLOSURE 

We have presented a robust and computationally efficient non-negative framework for fast 
bimolecular diffusive-reactive systems. We have rewritten the governing equations for the concen- 
tration of reactants and product in terms of invariants which are unaffected by the underlying 
reaction. Because of this algebraic transformation, instead of solving three coupled (nonlinear) 
diffusion-reaction equations, we need to solve only two uncoupled linear diffusion equations. This 
will considerably decrease the computational cost and also avoid numerical challenges due to non- 
linear terms. The main finding of the this paper is that a small violation of the non-negative 
constraint and maximum principle in the calculation of invariants can result in bigger errors in the 
prediction of the concentration of the product. 

Using several numerical experiments, it has been shown that the standard single-field formu- 
lation (which is also known as the Galerkin formulation) gives unphysical negative values for the 
concentration of the invariants, reactants, and the product. The formulation even gives unphysical 
negative value for the average of the concentration of the product, which is a useful measure in the 
prediction of the fate of transport of chemical species. It is also shown that the clipping procedure 
(which is considered as a variational crime) is not a viable fix to obtain the non-negative con- 
centrations, and does not give accurate results for diffusive-reactive systems. On the other hand, 
the proposed computational framework provided accurate results, and in all cases the framework 
provided non- negative values for the concentrations of the chemical species. The proposed com- 
putational framework can serve as a robust simulator for anisotropic heterogeneous geomodels. A 
next logical step is to extend the proposed computational framework to slow reactions, and to more 
complicated reactions. 
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Figure 1. This figure illustrates the violation of the comparison principle by the 
optimization-based solver. Note that the Galcrkin formulation also violates the comparison 
principle. The problem is a one-dimensional problem with homogeneous Dirichlet boundary 
conditions. The decay coefficient is taken as 2 x 10 4 , the diffusivity coefficient is taken as 
unity, and the domain is meshed using six linear finite elements. As evident from the figure, 
we have f 3 >z f 2 h f\ but the obtained numerical solution does not obey C3 y c 2 >r C\ . 
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Figure 2. Numerical /i-convergence study: A pictorial description of the boundary value 
problem. The computational domain is a rectangle with origin at x or i g in = (0,0). The 
length and width of the rectangular domain is L x — 2 and L y — 1. The volumetric sources 
for invariants F and G are prescribed. The Dirichlet boundary conditions c^(x) and cji,(x) 
are prescribed directly by evaluating the expressions given by equations (5.4a I and ( 5.4b| ) 
on the boundary. 




Figure 3. Numerical /i-convergence study: This figure shows the typical computational 
meshes (using three-node triangular and four-node quadrilateral elements) employed in the 
numerical convergence study. Hierarchical meshes are employed in the study. The meshes 
shown in the figure have 11 nodes along each side of the domain. Convergence analysis is 
performed using 11 x 11, 21 x 21, 41 x 41, and 81 x 81 meshes. 
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Figure 4. Numerical /i-convergence study: This figure shows the variation of error in the 
concentrations of invariants F and G with respect to mesh refinement. The left figure is for 
three- node triangular element, and the right figure is for four-node quadrilateral element. 
The error is calculated in both L2 norm and H 1 seminoma. Since the obtained values for the 
error are identical for both F and G, only the error for invariant F is shown in the figure. 




Figure 5. Numerical /i-convergence study: This figure shows the variation of error in the 
concentrations of A, B, and C with respect to mesh refinement. The left figure is for three- 
node triangular element, and the right figure is for four-node quadrilateral element. It is 
evident that the numerical solution obtained from the proposed computational framework 
converges with respect to mesh refinement. 
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FIGURE 6. Plume development from boundary in a reaction tank: A pictorial description 
of boundary value problem. The computational domain is a rectangle with reactants A and 
B are pumped on the left face of the domain, which are modeled as Dirichlet boundary con- 
ditions. On the remaining faces of the domain, no flux (i.e., Neumann) boundary condition 
is applied. Dirichlet boundary conditions for A and B are also indicated in the figure. 
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Figure 7. This figure shows the contours of the multi-mode stream function given by equa- 



tion (5.8), and corresponding vector field given by equation (5.9). It should be emphasized 
that the velocity should not be associated vMh the advection velocity, as advection is ne- 
glected in this paper. Herein, the velocity (given by equation ( |5.9| ) should be interpreted as 
the principal directions (or eigenvectors) of the diffusivity tensor (which is given by equation 



(5.7)). Since the streamlines do not intersect, the nature of the vector fields can provide 
insights on the plume of the reactants and the product. 
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Figure 8. Plume formation from boundary in a reaction tank: This figure shows the 
contours of the invariant F obtained using the Galerkin formulation (top figure) , the clipping 
procedure (middle figure), and the proposed computational framework (bottom figure). The 
concentration of the invariant F should be between and 1. The Galerkin formulation 
violates the non-negative constraint and the maximum principle. On the other hand, by 
explicitly enforcing the lower and upper bounds, the proposed computational framework 
satisfies the non-negative constraint and the maximum principle. However, it is noteworthy 
that the contours of the invariant F under the Galerkin formulation, the clipping procedure 
and the proposed computational framework all look similar. This is not the case with the 



contours of the product C (see Figure 12 ) 
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Figure 9. Plume formation from boundary in a reaction tank: This figure shows the 
contours of the invariant G obtained using the Galerkin formulation (top figure) , the clipping 
procedure (middle figure), and the proposed computational framework (bottom figure). The 
concentration of the invariant G should be between and 10. The Galerkin formulation 
violates the non-negative constraint and the maximum principle. On the other hand, by 
explicitly enforcing the lower and upper bounds, the proposed computational framework 
satisfies the non-negative constraint and the maximum principle. However, it is noteworthy 
that the contours of the invariant G under the Galerkin formulation, the clipping procedure 
and the proposed computational framework all look similar. This is not the case with the 



contours of the product C (see Figure 12). 
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Figure 10. Plume formation from boundary in a reaction tank: This figure shows the 
contours of the reactant A obtained using the Galerkin formulation (top figure), the clipping 
procedure (middle figure), and the proposed computational framework (bottom figure). The 
concentration of the reactant A should be less than or equal to 1. The Galerkin formulation 
violates the maximum principle. On the other hand, the proposed computational framework 
satisfies the upper bound due to the maximum principle. However, it is noteworthy that the 
contours of the reactant A under the Galerkin formulation, the clipping procedure and the 
proposed computational framework all look similar. This is not the case with the contours 



of the product C (see Figure 12 ) 
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Figure 11. Plume formation from boundary in a reaction tank: This figure shows the 
contours of the reactant B obtained using the Galerkin formulation (top figure), the clipping 
procedure (middle figure), and the proposed computational framework (bottom figure). The 
concentration of the reactant B should be less than or equal to 10. The Galerkin formulation 
violates the maximum principle. On the other hand, the proposed computational framework 
satisfies the upper bound due to the maximum principle. However, it is noteworthy that the 
contours of the reactant B under the Galerkin formulation, the clipping procedure and the 
proposed computational framework all look similar. This is not the case with the contours 



of the product C (see Figure 12). 
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Figure 12. Plume formation from boundary in a reaction tank: This figure shows the 
contours of the concentration of the product C obtained using the Galerkin formulation 
(top figure) , the clipping procedure (middle figure) , and the proposed computational frame- 
work (bottom figure). Mathematically, the concentration of the product C should be non- 
negative. There is no upper bound on the concentration of the product due to the maximum 
principle, as there is a source for the product C due to the underlying reaction. As one can 
see from the figure, the Galerkin formulation violates the non-negative constraint. On the 
other hand, the proposed computational framework satisfies the non-negative constraint 
for the concentration of the product C. More importantly, the contours of the product 
under the Galerkin formulation, the clipping procedure and the proposed computational 
framework differ considerably. Using this figure, one can draw an important conclusion of 
the paper: A small violation of the non-negative constraint and the maximum principle in 
the calculation of the invariants can result in significant errors in the concentration of the 
product for diffusive-reactive systems. „-, 



Correspondence to: Dr. Kalyana Babu Nakshatrala, Department of Civil and Environmental 
Engineering, Engineering Building #1, Room #135, University of Houston, Houston, Texas 77204- 
4003, USA. TEL: + l-713-743-4418 

E-mail address: knakshatrala@uh.edu 

Maruti Kumar Mudunuru, Graduate Student, Department of Civil and Environmental Engineer- 
ing, University of Houston, Houston, Texas 77204-4003, USA. 

Dr. Albert Valocchi, Department of Civil and Environmental Engineering, Newmark Civil Engi- 
neering Laboratory, University of Illinois at Urbana-Champaign, Urbana, Illinois 61801, USA. 



.38 



0.2 Q.4 0.6 0.6 1 




2 4 6 a 10 




0.00 0.07 0.14 0.22 0.29 0,36 




Figure 13. Plume formation from boundary in a reaction tank: This figure shows the 
contours of invariant F (top figure), invariant G (middle figure), and product C (bottom 
figure) under the Galerkin formulation. For the invariant F, the concentrations outside the 
interval [0, 1] are cut-off, and for the invariant G the concentrations outside [0, 10] are cut- 
off. For the product C, the concentrations below zero are cut-off. The lower bounds are due 
to the non-negative constraint, and the upper bounds are due to the maximum principle. 
The regions with cut-off concentration values are indicated in white color. It is evident that 
the Galerkin formulation violated the bounds in significant portions of the computational 
domain. 39 
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Figure 14. Plume formation from boundary in a reaction tank: This figure shows the 
variation of the concentration of the product at y = 0.45 (top figure) and x = 1.75 (bottom 
figure). The Galerkin formulation produces negative concentrations for the product that is 
extended over significant region, and the violation is not mere numerical noise. It should 
also be noted the concentration of the product under the clipping procedure is different from 
that of the proposed computational framework. 
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Figure 15. Plume formation from boundary in a reaction tank: This figure shows the 
variation of integrated concentration of the invariants F and G with respect to y along x 
(i.e., the variations of J cp (x, y) dy and Jcg(x, y) dy along x) for various meshes using the 
Galerkin formulation, and the proposed computational framework. The results obtained 
under the clipping procedure are almost the same as that of the Galerkin formulation and 
the proposed computational framework, and hence not shown in the figure. 
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Figure 16. Plume formation from boundary in a reaction tank: This figure shows the 
variation of integrated concentration of the product C with respect to y along x (i.e., the 
variation of J cc (x, y) dy along x) for various meshes using the Galerkin formulation, the 
clipping procedure, and the proposed computational framework. This figure clearly shows 
that the Galerkin formulation predicts unphysical negative values even for statistical average 
quantities. In particular, the Galerkin formulation predicts negative values for the cross- 
sectional total concentration of the product for significant portions along the x-direction 
even on fine computational grids. 
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Figure 17. Plume formation from boundary in a reaction tank: This figure shows the 
contours of the natural logarithm of the Lagrange multipliers corresponding to the lower 
bound constraints (top) and the upper bound constraints (bottom) for the invariant F. The 
proposed computational framework is employed, which solves a quadratic programming 
problems to obtain the concentration of invariant F. It should be noted that the Lagrange 
multipliers will always be non-negative, and hence the natural logarithm of the Lagrange 
multipliers will be (extended) real numbers. For visualization purposes, we have set the 
natural logarithm of zero to be —100. 
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Figure 18. Plume formation from boundary in a reaction tank: This figure shows the 
contours of the natural logarithm of the Lagrange multipliers corresponding to the lower 
bound constraints (top) and the upper bound constraints (bottom) for the invariant G. The 
proposed computational framework is employed, which solves a quadratic programming 
problems to obtain the concentration of invariant G. It should be noted that the Lagrange 
multipliers will always be non-negative, and hence the natural logarithm of the Lagrange 
multipliers will be (extended) real numbers. For visualization purposes, we have set the 
natural logarithm of zero to be —100. 
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Figure 19. Plume formation due to multiple stationary point sources: A pictorial de- 
scription of boundary value problem. The computational domain is a rectangle with ori- 
gin at x or igi n = (0,0). The length and width of the rectangular domain is L x = 2 and 
L y = 1. Continuous point sources of chemical species A and B are located at Xi = (0.2, 0.3), 
x 2 = (1.6,0.3), x 3 = (1.6,0.7), and X4 = (0.2,0.7) release reactants at a constant rate of 
f\ = 0.1, f\ = 0.05, fg = 0.1, and fg = 0.1 moles. The Dirichlet boundary conditions 
c^(x) , c'g(x), and c^(x) are prescribed to be zero on the boundary. 
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Figure 20. Plume formation due to multiple stationary point sources: This figure shows 
the steady-state contours of concentration of the invariant F under the Galerkin formulation 
(top figure), and the non- negative formulation (bottom figure). The minimum and maximum 
concentrations of F under the Galerkin formulation are respectively —0.0498 and 4.7100, 
and the corresponding values under the non-negative formulation are respectively and 
4.7088. The rotated heterogeneous diffusivity tensor is employed in the numerical simulation. 
The results are generated using 201 x 201 structured mesh using four-noded quadrilateral 
elements. The regions with negative concentration are indicated in white color. 
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Figure 21. Plume formation due to multiple stationary point sources: This figure shows 
the steady-state contours of concentration of the invariant G under the Galerkin formulation 
(top figure), and the non- negative formulation (bottom figure). The minimum and maximum 
concentrations of G under the Galerkin formulation are respectively —0.0158 and 0.9977, 
and the corresponding values under the non-negative formulation are respectively and 
0.9969. The rotated heterogeneous diffusivity tensor is employed in the numerical simulation. 
The results are generated using 201 x 201 structured mesh using four-noded quadrilateral 
elements. The regions with negative concentration are indicated in white color. 
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Figure 22. Plume formation due to multiple stationary point sources: This figure shows 
the steady- state contours of concentration of product C under the Galerkin formulation (top 
figure), and the non- negative formulation (bottom figure). The minimum and maximum 
concentrations of C under the Galerkin formulation are respectively —0.0996 and 0.2560, 
and the corresponding values under the non-negative formulation are respectively and 
0.2564. The rotated heterogeneous diffusivity tensor is employed in the numerical simulation. 
The results are generated using 201 x 201 structured mesh using four-noded quadrilateral 
elements. The regions with negative concentration are indicated in white color. 
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Figure 23. Plume formation due to multiple stationary point sources: This figure shows 
the variation of integrated steady-state concentration of the product C with respect to y 
along x (i.e., the variation of J cc(x, y) dy along x) for various meshes using the Galerkin 
formulation and the proposed computational framework. Three different structured meshes 
using four-node quadrilateral elements are employed, which are indicated as 101 x 101, 
161 x 161 and 201 x 201. The rotated heterogeneous diffusivity tensor is employed in the 
numerical simulation. 
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Figure 24. Plume formation due to multiple stationary point sources: This figure compares 
the CPU timing of the Galerkin formulation and the proposed computational framework. 
Note that XSeed denotes the number of nodes along x-direction, and same number of nodes 
are employed in the y-direction. In other words, a mesh with XSeed = YSeed = 201 will have 
over 40,000 total number of nodes. For each node, we have 5 unknonws - two invariants, two 
reactants, and the product. Note that we have employed MATLAB's interior-point-convex 
algorithm to solve the resulting quadratic programming problems. Each point is generated 
by taking the average of five simulations. 
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Figure 25. Diffusion and reaction of an initial slug: A pictorial description of the initial 
boundary value problem. The slug of dimensions d x = 2 , d y = 1 is located at x s i ug = (4, 2) 
in the rectangular domain of length L x = 10 and breadth L y = 5. 
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Figure 26. Diffusion and reaction of an initial slug: This figure shows the contours of the 
concentration of the product C obtained using the Galerkin formulation (top figure), and 
the proposed computational framework (bottom figure) at t = 0.05 s. The time step is 
taken as At — 0.05 s. As one can see from the figure, the Galerkin formulation violated 
the non- negative constraint. On the other hand, the proposed computational framework 
satisfies the non-negative constraint for the concentration of the product C. 
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Figure 27. Diffusion and reaction of an initial slug: This figure shows the variation of 
the concentration of the product at y = 2.5 at three time levels. The Galerkin formulation 
produces negative shoots for the concentration of the product, and the violation is not mere 
numerical noise. The proposed computational framework produced physically meaningful 
non-negative values for the concentration of the product. 
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Figure 28. Diffusion and reaction of an initial slug: This figure shows the number of 
iterations taken by the quadratic programming solver at each time level in obtaining the 
concentrations for the invariants. The time step is taken as 0.05. 
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